Mean Field

What You Will Learn

• The Mean Field Approximation (MFA) simplifies complex systems by replacing fluctuating variables with their average values.

• When applied to functions of many interacting variables, MFA assumes that correlations can be neglected:

\[ \langle f(x_1, x_2, x_3, \dots) \rangle \approx f(\langle x_1 \rangle, \langle x_2 \rangle, \langle x_3 \rangle, \dots) \]

• This approximation is foundational across many areas of physics, including:

  • Electronic structure theory (e.g., Hartree–Fock)

  • Liquid state theory (e.g., van der Waals)

  • Condensed matter physics (e.g., Ising and Heisenberg models)

MFA applied to the Ising Model

• In the Ising model, the mean-field approximation (MFA) decouples spin–spin interactions, reducing the system to a single macroscopic parameter: the average magnetization per spin, denoted by \( m \).

\[ \langle s_i s_j \rangle \approx \langle s_i \rangle \langle s_j \rangle = m^2 \]

\[ \langle H \rangle \approx -J \sum_{\langle ij \rangle} \langle s_i \rangle \langle s_j \rangle - h \sum_i \langle s_i \rangle = -\frac{1}{2} J q m^2 - h m \]

• The factor of \( \frac{1}{2} \) avoids double-counting spin pairs, and \( q \) is the coordination number—for instance, \( q = 4 \) for a 2D square lattice and \( q = 6 \) for a 3D cubic lattice.

Derivation of the Mean-Field Hamiltonian

• Consider the 2D Ising model where each spin \( s_i \) interacts with its \( q \) nearest neighbors:

\[ H_i = -J \sum_{\delta} s_i s_{i+\delta} - h s_i = -\left( J \sum_{\delta} s_{i+\delta} + h \right) s_i \]

• Define the local mean field at site \( i \) as:

\[ H_i^{\text{MF}} = J \sum_{\delta} s_{i+\delta} + h \]

• Under the mean-field approximation, replace neighboring spins with their average value:

\[ H_i^{\text{MF}} \approx J q m + h \]

• The total Hamiltonian then becomes:

\[ H \approx - (J q m + h) \sum_i s_i \]

• Taking the thermal average and accounting for double counting leads to the previously derived mean-field expression for \( \langle H \rangle \).

Self-Consistent Equation for Magnetization

• Partition Function: After applying MFA, the system becomes effectively non-interacting, and the total partition function factorizes as \( Z = z^N \), where \( z \) is the partition function of a single spin. For \( s_i = \pm 1 \):

\[ E(s_i = \pm 1) = \mp (J q m + h) \]

\[ z = e^{\beta (J q m + h)} + e^{-\beta (J q m + h)} \]

• The probability that a spin is in state \( s_i \) is:

\[ p(s_i) = \frac{e^{\beta s_i (J q m + h)}}{z} \]

• Average Magnetization: The mean value of a spin is then:

\[ m = \langle s_i \rangle = \sum_{s_i} s_i \, p(s_i) = \frac{e^{\beta (J q m + h)} - e^{-\beta (J q m + h)}}{e^{\beta (J q m + h)} + e^{-\beta (J q m + h)}} \]

\[ m = \tanh(\beta (J q m + h)) \]

• This self-consistent equation relates the average magnetization \( m \) to itself via the hyperbolic tangent function, capturing the thermal competition between spin alignment and disorder.

Free Energy of Mean-Field Ising Models

Defining the Macrostate \( M \)

• Consider a spin lattice of \( N \) sites, where each spin can point up or down. The total magnetization is given by:

\[ N = N_{+} + N_{-}, \qquad M = N_{+} - N_{-} \]

• Defining the magnetization per spin \( m = M/N \), the probabilities of spin-up and spin-down states are:

\[ p_{\pm} = \frac{1 \pm m}{2} \]

Entropy, Energy, and Free Energy

• Entropy is given by the Shannon expression using the probabilities of spin states:

\[ S(m) = -N k_B \left[ \left( \frac{1+m}{2} \right) \log \left( \frac{1+m}{2} \right) + \left( \frac{1-m}{2} \right) \log \left( \frac{1-m}{2} \right) \right] \]

• Energy in the mean-field approximation is the sum over spins in the presence of an average field:

\[ E(m) = -N \left[ \frac{Jq}{2} m^2 + h m \right] \]

• Free Energy is the usual Legendre transform:

\[ F(m) = E(m) - T S(m) \]

• Dimensionless Free Energy per Spin (useful for analysis and plotting):

\[ f(m) = \frac{F}{N k_B T} = - \left[ \frac{Jq}{2 k_B T} m^2 + \frac{h}{k_B T} m \right] + \left[ \left( \frac{1+m}{2} \right) \log \left( \frac{1+m}{2} \right) + \left( \frac{1-m}{2} \right) \log \left( \frac{1-m}{2} \right) \right] \]

Finding the Critical Temperature \( T_c \)

• The critical point is determined by the second derivative of the free energy:

\[ \frac{\partial^2 f}{\partial m^2} \bigg|_{m=0} = 1 - \frac{Jq}{k_B T} = 0 \]

\[ \Rightarrow \quad T_c = \frac{Jq}{k_B} \]

Magnetization as a Function of Temperature

• The equilibrium magnetization minimizes the free energy, leading to a self-consistency condition:

\[ m = \tanh \left( \frac{Jq}{k_B T} m + \frac{h}{k_B T} \right) \]

• Rearranged as:

\[ - \left( \frac{Jq}{k_B T} \right) m - \frac{h}{k_B T} + \tanh^{-1}(m) = 0 \]

• For zero external field (\( h = 0 \)), the equation simplifies to:

\[ m = \tanh \left( \frac{T_c}{T} m \right) \]

• For small \( m \), expand \( \tanh^{-1}(m) \approx m + \frac{1}{3} m^3 \), yielding:

\[ - \frac{T_c}{T} m + m + \frac{1}{3} m^3 = 0 \]

• Solving this, we find the magnetization near the critical point behaves as:

\[ m = \pm \sqrt{ \frac{3(T_c - T)}{T_c} } \]

Visualizing MF predictions

The \(h=0\) MFA case

The equation can be solved in a self-consistent manner or graphically by finding intersection between:

• \(m =tanh(x)\)

• \(x = \frac{Jqm}{k_BT}\)

When the slope is equal to one it provides a dividing line between two behaviours.

\[k_B T_c =qJ\]

\[m = tanh \Big(\frac{Tc}{T} m \Big)\]

MFA shows phase transitio for 1D Ising model at finite \(T=T_c\)!

import ipywidgets as widgets
from ipywidgets import interactive, interact
import matplotlib.pyplot as plt
import numpy as np

# Constants
J = 1  # Interaction strength
k = 1  # Boltzmann constant (set to 1 for simplicity)

def entropy(m):
    # Handle log of zero by adding a small number to the argument
    small_number = 1e-10
    return -(1+m)/2 * np.log((1+m)/2 + small_number) - (1-m)/2 * np.log((1-m)/2 + small_number)

def free_energy(m, T):
    # Calculate the free energy for given m and T
    return -0.5 * J * m**2 - T * entropy(m)

def plot_free_energy(T):
    # Magnetization range
    m_values = np.linspace(-1, 1, 100)
    F_values = [free_energy(m, T) for m in m_values]
    
    # Plotting
    plt.figure(figsize=(8, 5))
    plt.plot(m_values, F_values, label=f'T = {T}')
    plt.xlabel('Magnetization (m)')
    plt.ylabel('Free Energy per Spin (F)')
    plt.title('Mean-field Free Energy vs. Magnetization')
    plt.legend()
    plt.grid(True)
    plt.show()

#interactive(plot_free_energy, T=(0.1, 2, 0.1 ))

def mfa_ising_Tc(T=1, Tc=1):

    x  = np.linspace(-3,3,1000)
    
    f = lambda x: (T/Tc)*x
    m = lambda x: np.tanh(x)
    
    plt.plot(x,m(x), lw=3, alpha=0.9, color='green')
    
    plt.plot(x,f(x),'--',color='black')
    idx = np.argwhere(np.diff(np.sign(m(x) - f(x))))
    plt.plot(x[idx], f(x)[idx], 'ro')
    
    plt.legend(['m=tanh(x)', 'x'])
    
    plt.ylim(-2,2)
    plt.grid('True')
    plt.xlabel('m',fontsize=16)
    plt.ylabel(r'$tanh (\frac{Tc}{T} m )$')
    plt.show()

#interact(mfa_ising_Tc, T=(0.1,5)) 

def mfa_ising_h_vs_m(T=1):

    Tc = 1
    x  = np.linspace(-1,1,200)
    
    h = T*(np.arctanh(x) - (Tc/T)*x)
    
    plt.plot(x, h, lw=3, alpha=0.9, color='green')

    plt.plot(x, np.zeros_like(x), lw=1, color='black')
    plt.plot(np.zeros_like(x), x, lw=1, color='black')
     
    plt.grid(True)

    plt.ylabel('m',fontsize=16)
    plt.xlabel('h',fontsize=16)
    plt.ylim([-1,1])
    plt.xlim([-1,1])
    plt.show()

#interactive(mfa_ising_h_vs_m, T=(0.1,5)) 

import plotly.graph_objects as go
import numpy as np
from scipy.optimize import fsolve
from scipy.optimize import root_scalar  # Importing root_scalar

def compute_xcross(T, h_over_T):
    
    def f(M):
        Tc = 1.0  # Normalized temperature
        h = h_over_T * T  # Compute actual h from h/T and T
        return M - np.tanh((M + h) / (T / Tc))
    
    # Using symmetric interval for root finding to allow negative solutions
    interval = [-2.0, 2.0]

    # Check if a root is likely within the given range
    if f(interval[0]) * f(interval[1]) > 0:
        return 0.0  # If no root likely, return 0
    
    # Find the root using bisection method within the specified range
    result = root_scalar(f, bracket=interval, method='bisect')
    
    return result.root 

# Define ranges for h/T  and T/Tc (from 0.2 to 1.5)
h_over_T_values = np.linspace(-0.2, 0.2, 101)
T_over_Tc_values = np.linspace(0.1, 1.6, 101)

# Create a meshgrid for T/Tc and h/T
H_over_T, T_over_Tc = np.meshgrid(h_over_T_values, T_over_Tc_values)

# Initialize an array for magnetizations
Magnetizations = np.zeros_like(H_over_T)

# Compute M for each (h/T, T/Tc) pair in the meshgrid
for i, T in enumerate(T_over_Tc_values):
    for j, h in enumerate(h_over_T_values):
        
        Magnetizations[i, j] = compute_xcross(T, h)


# Create a 3D surface plot using Plotly
fig = go.Figure(data=[go.Surface(z=Magnetizations, x=H_over_T, y=T_over_Tc, 
                                 colorscale='RdBu', 
                                 contours={
                                     #'z': {'show': True, 'start': -1.0, 'end': 1.0, 'size': 0.1, 'color':'orange'},
                                     'x': {'show': True, 'start': -1, 'end': 1, 'size': 0.05, 'color':'black'},
                                     'y': {'show': True, 'start': 0.5, 'end': 1.5, 'size': 0.1, 'color':'black'}
                                 })])

# Update plot layout
fig.update_layout(
    title='3D Plot of Magnetization M vs normalized temperature and field; T/Tc and h/T',
    scene=dict(
        xaxis_title='h/T',
        yaxis_title='T/Tc',
        zaxis_title='M',
        aspectmode='manual',
        yaxis_autorange='reversed',  # Reverse the Y-axis
        aspectratio=dict(x=2, y=2, z=0.75),  # Adjust these values as needed
        camera=dict(
            eye=dict(x=-4, y=-2, z=2),  # Adjust camera "eye" position for better view
            up=dict(x=0, y=0, z=1)  # Ensures Z is up
        )
    ),
    autosize=False,
    width=900,
    height=900
)

# Display the figure with configuration options
fig.show()

plt.plot(Magnetizations[:, len(h_over_T_values)//2-1], color='blue')
plt.plot(Magnetizations[:, len(h_over_T_values)//2+1], color='blue')

[<matplotlib.lines.Line2D at 0x7f90a447f580>]